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AN  EVALUATION  OF  THE  METHOD  OF  DETERMINING 
PARALLAX  FROM  MEASURED  PHASE  DIFFERENCES 


INTRODUCTION 

The  purpose  of  the  report  is  to  describe  an  evaluation  of  the  practicality  of  deter- 
mining parallax  by  means  of  detecting  phase  differences  extracted  from  corresponding 
digital  signals.  The  plan  was  to  perform  a preliminary  evaluation  by  registering  an  image  to 
itself  using  measured  phase  differences  and  if  the  test  results  proved  promising,  the  plan  was 
to  enlarge  the  scope  of  the  tests  to  include  stereo  images.  Stereo  tests  were  not  performed 
because  the  initial  results  demonstrated  that  conventional  autocorrelation  methods  were 
more  accurate  and  more  practical  than  measuring  parallax  by  comparing  phase  differences. 

The  term  autocorrelation  in  this  report  pertains  to  the  process  of  establishing  corre- 
spondence between  identical  image  pairs.  Conventional  digital  autocorrelation  processes 
involve  determining  maximum  or  minimum  points  of  discrete  functions.  The  functions  are 
generated  by  using  any  one  of  several  measures  of  simularity,  e.g.  the  linear  correlation  co- 
efficient described  in  the  Appendix.  The  term  crosscorrelation  pertains  to  the  process  of 
establishing  correspondence  between  dissimilar  images,  e.g.  stereo  images.  Conventional 
digital  crosscorrelation  processes  are  identical  to  conventional  autocorrelation  processes 
except  that  input  data  is  extracted  from  corresponding  scenes  rather  than  from  identical 
scenes. 


GEOMETRIC  DESCRIPTION 

In  the  proposed  evaluation,  corresponding  lines  of  image  data  are  extracted  from 
stereo  exposures,  and  specific  points  on  the  lines  are  brought  into  corresporxlence  by  signal 
matching  in  frequency  space. 

Epipolar  Geometry.  Consider  two  overlapping  frame  camera  exposures.  An  epipolar 
line  pair  results  from  the  intersection  of  the  associated  epipolar  plane  and  the  two  focal 
planes.  An  epipolar  plane  is  any  member  of  the  family  of  planes  defined  by  the  straight  line 
connecting  the  two  exposure  positions . The  corresponding  epipolar  section  is  the  planar 
space  curve  defined  by  the  epipolar  plane  and  the  surface  of  the  earth.  Then,  except  for 
points  obscured  by  terrain  heights,  a one-to-one  correspondence  exists  between  the  image 
points  on  epipolar  line  pairs  and  the  associated  object  space  point  found  on  the  epipolar 
section  (See  figure  1). 

Suppose  that  Pi  is  selected  on  the  first  image  and  that  the  epipolar  lines  are  con- 
structed, then  the  corresponding  image  point  Pi  is  somewhere  on  the  second  epipolar 
line.  Thus,  if  the  epipolar  lines  can  be  constructed,  only  a one-dimensional  search  need  be 
conducted  for  Pi . Epipolar  lines  can  be  determined  if  the  relative  orientation  of  the  two 
exposures  is  known. 


5 


CSI  O C>J 


AX 

L 


POINT  DEFINED  ON  IMAGE  #1 
CORRESPONDING  POINT  ON  IMAGE  12 
ASSUMED  MATCH  POINT 
P2  - Pf  : X-PARALLAX 
LENGTH  OF  SCAN 


FIGURE  2.  EP I POLAR  SCANS. 
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Fourier  Representation.  Suppose/gt  (i):  i = - (L-a  )/j , L/j)}  is  a set  of  L equally 
?paced  density  values  centered  about  Pi  on  tbe  first  epipolar  line.  Suppose  also  that 
(g»  (i);  i = - (L-a  )/a , L/a^is  a set  of  L equally  spaceo  jensity  values  centered  about  (an 
estimate  of  Pa  on  the  corresponding  epipolar  line.  The  sit  uation  is  described  in  figure  2 
below,  xhe  fourier  expansion  of  the  two  signals  is  as  follows: 


gj(i)  * .^  + 


82(1)  • 


yk"  l/«K  + ' 

«k  • f ' Tan 

(aj^,  b|^):  Harmonic  coefficients  from 
discrete  Fourier  transform 


The  Fourier  representation  requires  that  the  signals  be  periodic,  they  are  not;  and  the 
match  technique  implies  that  the  two  signals  are  simply  translated  with  respect  to  one 
another,  they  are  not.  These  assumptions,  and  others,  are  discussed  in  the  section  on 
evaluation. 


Phase  Shift  and  Parallax.  From  figure  1,  the  two  curves  can  be  brought  into 
registration  by  translating  the  origin  of  gj  (i)  by  A X = Pj  - P2  or  by  letting  i = i ' + 
AX. 


1£  + S Yk  Sin{  2Kir  {i'+  AX)  + 6^^ 


L 


«K2! 
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then 


If  gi(0  ' 

6K:  • * iKz 

L 

for  each  harmonic  K. 


The  paraJJax  value  aX  can  be  estimated  from  each  of  the  L harmonics.  A weighted 
estimate  of  hX  will  be  calculated  from  Q of  the  harmonics. 


AX 


Wq  AXq 

Wq  = 1 


The  weights  will  be  defined  in  the  next  section. 


NUMERICAL  EXPERIMENTS 

The  digitued  image  used  for  this  experiment  is  stored  in  the  DlAL^Digital 
Image  Analysis  Laboratory)  system,  and  it  is  described  in  a previous  ETL  repoit.  The 
digital  image  is  a subimage  taken  from  an  exposure  of  rugged  terrain.  The  essential 
features  of  the  digitized  model  are 

Exposure  Data 

Camera  Type.  Vertical  Frame 
Focal  Length;  6 inches  (15.2  cm) 

Scale  1:48000 


*UwT«nce  A.  Gtmbino  ind  Bryce  L Schlock,  “An  Experimenttl  Digitil  Inlernctive  Ficiity”,  Compoler,  Vol.  10. 
No.  8.  Auiueet  1977,  pp.  22-28. 

^Michael  CiomMe,  Stereo  Anilyris  of  A SpecifK  Digital  Modei  Sampled  From  Aerial  Imagery,  U.S.  Army  Engineer 
Topographic  Laboratoriet,  Fort  Beivoir,  Va.,  ETL'0072,  September  1976  AD-A033  567. 
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Digital  Data 


Image  Dimension;  2048  x 2048 

Grayshades  ; 256 

Spot  Diameter  : 34.5  yifl  (micrometers) 

Spot  Spacing  24  pm 


The  experiment  was  one  of  autocorrelation;  therefore,  there  was  no  need  to 
extract  epipolar  line  data.  Approximately  25  lines  of  data  were  used  in  the  tests.  The 
lines  were  chosen  to  reflect  a large  range  in  signal  power.  Essentially,  a line  of  (L  + 2S) 
pixel  values  was  extracted  from  the  digital  image  on  disc.  The  center  of  the  line  was 
designated  as  the  test  point  Pi ; the  image  point  Pi  wac  characterized  by  a window  of  L 
pixels  about  K.  The  corresponding  point  Pi  was  defined  to  be  at  location  (Pi  + s), 
where  s varied  from  zero  to  S.  The  image  point  Pi  was  also  characterized  by  a window 
of  L pixels  about  it.  The  match  technique  was  applied  to  the  two  sets  of  pixels  to 
determine  how  accurately  the  induced  shift  of  + s could  be  recovered. 


Weiriits.  Two  weighting  functions  were  tested  in  this  experiment.  The  first 
was  derived  from  least  squares,  and  the  second  was  an  impirical  one  that  was  simpler 
to  compute.  It  can  be  shown  that  if  the  harmonic  coefficients  are  computed  by  least 
squares,  then  the  standard  error  of  the  parameter  6)^  is 


a6)(  » 

* \l4 

iTq’.  Standard  error  of  the  noise  associated 
with  the  g(iK 

N : Number  of  data  points,  L in  this  case. 


If  the  noise  is  the  same  for  both  scans,  then  it  can  be  also  show  that  the  standard 
error  of  the  parallax  is 


oAXk 


2 JT  [_  00 

W Y'' 
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If  the  parallax  value  AX  is  determined  by  leaat  iquaret  from  Q of  the  AXj^  values  and 
if  the  AXj^  are  weighted  according  to  oAXj^  above,  then  the  individual  weights  are 


where  q is  the  specific 
harmonic  number. 


The  second  weight  function  was  defined  to  be 


2 


0 


Window  Size.  The  window  size  L was  allowed  to  vary,  as  was  the  input  shift 
s. 


! « 16,  32,  64,  128,  256  pixels 
s»0,  + !,-••  2,  . . pixels 


Note  that  the  window  size  was  increased  in  length  and  the  pixel  spacing  remained  the 
same.  A fast  Fourier  transform  routine  from  the  DIAL  system  was  used  to  compute 
the  required  transform  data. 
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EVALUATION 


If  two  sinusoids  are  out  of  phase  by  period,  then  the  cross-correlation  is 
zero,  and  the  pull-in  range  can  be  represented  by  the  width  of  the  correlation  func- 
tion. Table  1 represents  PR,  the  pull-in  range  for  the  image  used  in  the  tests. 


TABLE  1.  PULL  IN  RANGE 


Ip/mm  PR  Pixels 


1 

+ 250 

+ 10.0 

2 

+ 125 

+ 5.0 

5 

+ 50 

+ 2.0 

7 

+ 35 

+ 1.5 

10 

+ 25 

+ 1.0 

20 

+ 13, 

+ 0.5 

Generally,  the  lower  frequencies  have  more  signal  power  associated  with  them  and, 
therefore,  carry  more  weight  in  the  shift  calculations.  Because  of  their  lower  power, 
the  higher  frequencies  are  more  adversely  affected  by  noise  than  are  the  lower  fre- 
quencies. The  fact  the  signal  and  its  corresponding  shifted  signal  are  not  periodic 
introduces  a distortion  into  the  calculations.  If,  in  the  stereo  case,  unshaped  data 
are  used,  then  an  additional  distortion  will  be  introduced. 


Fourier  Model.  The  fact  that  image  signals  are  not  periodic  caused  most  of 
the  inconsistency  in  the  method.  The  consistency  of  the  computed  shifts  AXj^ 
turned  out  to  be  highly  dependent  upon  the  consistency  of  the  signal  power  SP,  a« 
the  window  was  displaced  + s in  the  pixel  direction.  Signal  power  is  defined  here  as 
the  mean  squared  variation  of  the  AC  terms,  or  simply  as  the  variance  of  the  signal 

y (g  - g)^ 

SP  = -= 

L - 1 
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If  there  was  little  or  no  change  in  SP  as  the  induced  shift  s varied  over  its  range, 
then  the  match  results  were  accurate  and  precise.  If  SP  changed  by  10  percent  or 
more,  the  computed  varied  significantly  in  precision,  and  in  many  cases 

appeared  to  be  erratic. 


Stereo.  Consider  a ramp-like  object  imaged  on  two  overlapping  vertical  frame 
images.  If  a represents  the  tilt  of  the  ramp  with  respect  to  the  datum  and  B/H  repre- 
sents the  imaging  base-height  ratio,  then  the  entries  in  table  2 represent  scale  differ- 
ences between  the  two  images  when  the  object  is  located  in  the  middle  of  the  model. 
For  example,  if  a = 20®  and  B/H  = 1,  then  the  ramp  image  length  on  one  exposure  is 
1.44  times  the  i..  age  length  on  the  second  exposure.  It  is  assumed  here  that  the  ramp 
is  ahgned  in  the  epipolar  direction.  Mere  examples  of  scale  distortion  are  given  by 
Crombie  and  Gambino.^ 


TABLE  2.  IMAGE  SCALE  CHANGE  AT  MIDDLE  OF  MODEL 

a 


B/H 

0° 

5° 

10° 

15° 

O 

o 

CSI 

25° 

30° 

0.4 

1.00 

1.04 

1.07 

1.11 

1.16 

1.21 

1.26 

0.6 

1.00 

1.05 

i.n 

1.17 

1.25 

1.33 

1.42 

0.8 

1.00 

1.07 

1.15 

1.24 

1.34 

1.46 

1.60 

1.0 

1.00 

1.09 

1.19 

1.31 

1.44 

1.61 

1.81 

As  shown  in  table  2,  unless  one  image  is  stretched  to  accommodate  the  scale  change,  a 
large  amount  of  signal  difference  for  comparable  scan  lengtlis  will  occur,  especially 
for  large  values  of  a and  B/H. 

Numerical  Results.  The  results  given  in  table  3 pertain  to  Wq(i ),  the  first 
weighting  method.  This  procedure  turned  out  to  be  only  slightly  better  than  Wqfi). 
The  results  are  presented  as  a Ametion  of  Q and  L,  where  Q is  the  number  of  har- 
monics used  in  the  estimate  and  L is  the  window  length.  It  should  be  noted  that  it  is 
not  entirely  correct  to  compare  results  over  L.  (See  table  3). 


^Michael  A.  Crombie  and  Uwrence  A.  Gambino.  “Digital  Stereo  Photogrammetiy”,  U.S.  Army  Engineer  Topo- 
graphic Laboratories,  Fort  Belvoir,  Virginia,  Presented  at  the  Congress  of  the  International  Federation  of  Surveyo-  s 
Commission  V,  Stockholm,  Sweden,  June  1977. 
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TABLE  3.  LINE  PAIRS  PER  MILLIMETER 


H 


L 

1 

3 

5 

7 

16 

2.8 

8.3 

13.9 

19.5 

32 

1.3 

4.0 

6.7 

9.4 

64 

0.7 

2.0 

3.3 

4.6 

128 

0.3 

1.0 

1.7 

2.3 

256 

0.2 

0.5 

0.8 

1.1 

The  entries  in  table  3 are  line  pairs  per  millimeter  (1  p/mm).  They  were  calculated  by 
converting  the  window  lengths  into  millimeters  and  then  computing  the  number  of 
complete  cycles  per  millimeter  as  a function  of  the  harmonic  number  H.  For 
example,  the  third  harmonic  pertains  to  4 Ip  per  mm  when  L = 32  and  to  1 Ip  per 
mm  when  L = 128.  Thus,  when  Q = 7 and  L = 32,  resolution  infonnation  from  1.3  to 
9.4  Ip  per  mm  is  used  to  measure  parallax;  whereas,  when  Q = 7 and  L = 128,  the  reso- 
lution infonnation  ranges  from  0.3  to  2.3  Ip  per  mm. 

Tables  4 through  8 represent  sample  standard  errors  of  mismatch.  The  tabular 
parameters  are  defined  as 


S:  True  mismatch  in  pixels. 

Q;  Number  of  harmonics  used  in  the  weighted  average.  The  first  weight 
method  (modified  least  squares)  was  used.  For  example,  if  0 = 5, 
then  the  shift  was  estimated  from  the  first  5 harmonics. 

L:  Length  of  the  window  in  pixels. 

o.  Estimated  standard  error  with  approximately  50  degrees  of  freedom. 
The  standard  error  is  expressed  in  pixels. 
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TABLE  4.  O’  For  L * 16 


Q 


s 

1 

3 

5 

7 

+ 1 

0.76 

0.49 

0.47 

0.58 

s 

1 

TABLE  5.  0 For  L » 

Q 

3 

32 

5 

7 

+ 1 

0.59 

0.38 

0.32 

0.35 

+ 2 

1.22 

0.80 

0.73 

0.81 

s 

1 

TABLE  6.  «•  For  L = 64 

Q 

3 

5 

7 

+ 1 

0.50 

0.40 

0.36 

0.33 

+ 2 

0.91 

0.76 

0.71 

0.65 

+ 3 

1.33 

1.08 

1.07 

0.99 

4 4 

1.84 

1.40 

1.50 

1,66 

15 


TABLE  7.  #"  For  L » 128 


Q 


s 

1 

3 

5 

7 

+ 1_ 

0.53 

0.27 

0.17 

0.18 

+ 2 

1.04 

0.57 

0.37 

0.39 

+ 3 

1.68 

0.89 

0. 56 

0.60 

+ 4 

2.03 

1.19 

0.74 

0.78 

+ 5 

2.40 

1.47 

0.92 

1.04 

+ 6 

2.43 

1.73 

1.08 

1.20 

+ 7 

2.58 

1.92 

1.23 

1.64 

S 

1 

TABLE  8.  *r  For  L - 

Q 

3 

256 

5 

7 

+ 1 

0.64 

0.38 

0.26 

0.15 

+ 2 

1.41 

0.73 

0.52 

0.30 

± 3 

2.14 

1.07 

0.81 

0.46 

+ 4 

2.56 

1.25 

1.13 

0.62 

+ 5 

2.59 

1.32 

1.47 

0.83 

+ 6 

2.61 

1.59 

1.54 

1.00 

4 7 

2.71 

1.78 

1.76 

1.16 

16 


The  quantity  a/lsl  remains  nearly  constant  as  S varies  and  as  Q is  held  con- 
stant. This  observation  is  especially  true  for  L ^ 128  and  is  reasonably  true  for  L = 
256.  Tables  9 through  13  demonstrate  this  observation.  The  value  at  the  bottom  of 
each  column  is  the  geometric  mean  of  the  tabular  values. 

TABLE  9.  tf/lsl  For  L » 16 

Q 

S 1 3 5 7 

+ 1 0.76  0.49  0.47  0.58 


TABLE  10.  (T/lsl  For  L = 32 
Q 


s 

1 

3 

5 

7 

+ 1 

0.59 

0.38 

0.32 

0.35 

+ 2 

0.61 

0.40 

0.32 

0.41 

0.60 

0.39 

0.32 

0.38 

1 


TABLE  11.  <r/|s(For  L « 64 


Q 


s 

1 

3 

5 

7 

+ 1 

0.50 

0.40 

0.36 

0.33 

+ 2 

0.46 

0.38 

0.36 

0.33 

+ 3 

0.44 

0.36 

0.36 

0.33 

+ 4 

0.46 

0.35 

0.38 

0.42 

0.47 

0.37 

TABLE  12.  <r/<s|ForL 

Q 

0.37 

* 128 

0.35 

S 

1 

3 

5 

7 

■ 

+ 1 

0.53 

0.27 

0.17 

0.18 

+ 2 

0.52 

0.29 

0.19 

0.20 

+ 3 

0.56 

0.30 

0.19 

0.20 

+ 4 

0.51 

0.30 

0.19 

0.20 

+ 5 

0.48 

0.29 

0.18 

0.21 

+ 6 

0.41 

0.29 

0.18 

0.20 

+ 7 

0.37 

0.27 

0.18 

0.23 

0.49 

0.29 

0.18 

0.20 
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TABLE  13.  tf/lsl  For  L - 256 


Q 


s 

1 

3 

5 

7 

+ 1 

0.64 

0.38 

0.26 

0.15 

1 2 

0.71 

0.37 

0.26 

0.15 

+ 3 

0.71 

0.36 

0.27 

0.15 

+ 4 

0.64 

0.31 

0.28 

0.16 

1 5 

0.52 

0.26 

0.29 

0.17 

+ 6 

0.44 

0.27 

0.26 

0.17 

+ 7 

0.39 

0.26 

0.25 

0.17 

0.59 

0.32 

0.27 

0.16 

Coniider  the  average  values  of  a/lsl  in  table  14  and  the  Ip  per  mm  relation- 
ship to  specific  harmonic  values  listed  in  table  3.  A comparison  of  the  two  tables 
reveals  that  useful  parallax  information  can  be  derived  from  image  resolution  content 
ranging  from  approximately  1 to  about  10  Ip  per  mm. 

TABLE  14. 

AVERAGE  VALUES  OF  (r/{si 

L 

1 

Q 

3 

5 

7 

16 

0.76 

0.49 

0.47 

0.58 

32 

0.60 

0.39 

0.32 

0.38 

64 

0.47 

0.37 

0.37 

0.35 

128 

0.49 

0.29 

0.18 

0.20 

256 

0.59 

0.32 

0.27 

0.16 
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If  subpixel  precision  is  required  and  if  the  expected  parallax  S is  up  to  7 
pixels  or  larger,  then  an  averaging  of  parallax  estimates  at  and  around  the  point  in 
question  must  be  performed.  This  appears  to  be  a valid  approach  owing  to  the  ob- 
served variability  of  parallax  estimates  when  a line  is  shifted  in  the  pixel  direction. 
This  means  that  parallax  estimates  from  adjacent  lines  should  tend  to  be  independent. 
The  averaged  parallax  then  will  be  the  average  shift  over  the  area  defined  by  the  pixels 
used  in  the  computation.  Suppose  M lines  are  averaged,  the  standard  error  of  the 
parallax  is  required  to  be  f * (Pixel  Spacing)  where  0 ^ ^ ^ 1 , then  M must 
satisfy  the  following  relation. 


M 

and  suppose 

E: 

S: 

For  example,  suppose  pixel,  L = 64,  S = + 4 and  Q = 7,  then  E = 0.35  and  M 
> 31  lines.  If  in  the  same  example  S =^7,  then  M > 86  lines. 


is  the  appropriate  entry 
from  Table  14 

is  the  expected  shift 


The  results  of  this  report  pertain  to  autocorrelation.  It  has  been  shown  that  a 
large  number  of  calculations  need  to  be  performed  to  calculate  the  induced  shifts 
precisely.  Even  with  identical  images,  an  iterative  approach  must  be  performed  to 
reduce  the  parallax  error  to  zero.  It  should  be  noted  that  the  computations  can  be 
performed  efficiently  in  that  algorithms  can  be  developed  where  only  those  harmonics 
of  interest  need  be  computed.  It  should  also  be  noted  that  if  the  same  procedure 
were  to  be  applied  to  stereo  data,  then  the  following  error  sources  will  increase  the 
output  error: 

1 . Comparing  phase  angles  from  different  epipolar  lines.. 

2.  Comparing  phase  angles  when  one  line  is  stretched  with  respect  to  the  other. 

The  first  error  source  derives  from  erroneous  interior  and  exterior  orientation  data. 
The  second  error  source  occurs  whenever  the  corresponding  lines  are  not  shaped  accu- 
rately. When  the  parallax  is  significant,  the  second  error  source  has  the  most  effect  on 
the  output  error. 
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The  evaluation  described  in  this  report  must  be  performed  iteratively,  possibly 
with  diminishing  window  lengths  L,  to  achieve  zero  error  even  in  autocorrelation.  If 
accurate  and  precise  parallax  estimates  for  stereo  are  to  be  computed  in  relatively  steep 
terrain,  then  pixel  shaping  must  be  performed,  with  enough  care  taken  in  the  scan  pro- 
cedure to  ensure  that  corresponding  epipolar  lines  are  measured.  If  enough  care  is 
taken  to  produce  corresponding  epipolar  lines  and  if  the  lines  are  shaped  to  reflect 
expected  parallax,  then  conventional  correlation  methods  are  superior  to  frequency 
space  matching.  These  are  two  reasons  for  this.  First,  correct  shifts  can  be  computed 
from  a single  correlation  function  in  autocorrelation.  Second,  the  algorithm  (for 
stereo  matching  and  for  autocorrelation)  can  be  organized  so  that  every  point  on  corre- 
sponding epipolar  lines  can  be  matched  concurrently  with  great  efficiency.  The  appen- 
dix contains  the  computer  operations  associated  with  this  correlation  procedure.  A 
conventional  correlation  procedure  using  methods  very  similar  to  those  in  the  appendix 
has  been  applied  to  the  AS-1 1 B-X  device  that  was  developed  by  Bendix  Corporation 
for  the  Defense  Mapping  Agency.  The  great  economy  in  computer  operations  is 
derived  from  the  existence  of  significant  quantities  of  identical  calculations  as  the 
correlation  process  ranges  over  neighboring  points.  There  is  not  a similar  economy  in 
the  method  described  in  this  report.  The  utility  of  matching  all  possible  points  along 
the  corresponding  epipolar  pair  is  realized  when,  as  in  practice,  part  of  the  line  cannot 
be  matched  because  of  poor  image  detail  or  because  of  obscurations.  The  abundance 
of  matched  points  can  be  used  to  interpolate  for  missing  points  and  to  provide  a 
generally  smoother  set  of  output  points. 


CONCLUSIONS 


1 . Most  of  the  useful  parallax  information  can  be  derived  from  image  resolution  con- 
tent ranging  from  1 to  10  Ip  per  mm  for  the  aerial  image  used  in  this  report. 

2.  Measuring  parallax  by  comparing  phase  differences  between  corresponding  epi- 
polar lines  is  greatly  dependent  upon  stability  in  signal  power  as  the  matching  line  is 
shifted  in  the  pixel  direction. 

3.  If  signal  power  changed  by  10  percent  or  more,  the  computed  parallax  estimate 
varied  significantly  in  precision  and  in  many  cases  appeared  to  be  erratic. 

4.  Measuring  parallax  by  comparing  phase  differences  is  not  as  precise  as  conven- 
tional correlation  methods  in  autocorrelation.  It  is  expected  to  be  even  less  precise 
in  the  stereo  mode. 

5.  Measuring  parallax  by  comparing  phase  differences  does  not  offer  computational 
advantages  over  conventional  correlation  methods. 


^F.  Scarano  and  G.  Biumm  “A  Digital  Elevation  Data  Collection  Syitem",  PhotogrammetTic  Engineer  and  Remote 
Senaing.  Vol  XLII,  No.  4,  April  1976. 
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APPENDIX. 


Number  of  Computer  Operations  Needed  For  Vector  Correlation 


Given  two  N-dimensional  corresponding  epipolar  lines  of  gray  shades  X and  Y, 
suppose  M adjacent  points  of  Y are  to  be  matched  to  M corresponding  points  of  X. 
Assume  too  that  X has  been  shaped  approximately  so  that  the  M estimated  match 
points  of  X are  also  adjacent  points  or  nearly  so.  The  correspondence  will  be  refined 
by  correlation  methods  where  for  each  of  the  M matches,  P correlation  values  will  be 
calculated  using  a window  of  length  W.  If  P and  W are  odd  and  if  all  possible  points 
arc  matched,  then  M = N+  2-W-P. 

Essentially,  M discrete  correlation  functions  of  dimension  P will  be  developed 
where  the  central  value  of  each  pertains  to  the  approximate  match-point  location.  A 
shift  S for  each  of  the  M approximations  will  be  estimated  by  determining  the  loca- 
tion of  the  peak  correlation  value.  The  window  length  W must  be  large  enough  to 
produce  precise  correlation  values,  and  the  number  of  correlation  values  P must  be 
large  enough  so  that  the  extent  of  the  correlation  function  covers  the  actual  match 
location.  If  the  largest  possible  shift  in  absolute  value  is  SMAX,  then  (P-1)/-)  > 
SMAX.  The  purpose  of  this  appendix  is  to  calculate  the  number  of  computer  opera- 
tions necessary  to  calculate  the  M correlation  functions. 


The  measure  of  similarity  used  is  the  following  linear  correlation  coefficient; 


Where 

1 « 1,  M 

j = 1.  P 
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Any  one  of  the  required  M * P correlation  values  can  be  computed  from  the  following 
five  basic  quantities: 


and 


I 

k=1 


Y * Y 

^ijk  ’ijk 


A straightforward  count  of  the  necessary  number  of  computer  operations  is  givei' 
first.  No  regard  for  computational  overlap  is  taken  into  account.  However,  the  next 
calculation  pertains  to  the  situation  where  all  possible  computational  overlap  is  con 
sidered. 


Straightforward  Count.  The  subscripts  on  X and  Y are  dropped  for  con- 
venience. The  following  counts  of  the  number  of  add  and  multiply  operations  pertain 
to  the  five  basic  quantities  given  above. 


X : 

W adds  ^ 

X^: 

1 

W adds  J 

W multipl  ies  j 

1 

Y : 

> 

W adds 

r 

, 

W adds 

W multipl iesy 

I 

-XY: 

W adds  \ 

W multiplies  , 

* M * P 


★ M 
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The  total  number  of  adds  and  multiplies  for  the  five  components  is 


I 

! 


= M * W * (3  * P + 2) 

* M * W * (2  * P + 1) 

The  five  components  are  combined  in  the  following  way  to  produce  the  correlation 
value  r; 


wX  XY  -Zx  Z 


VwZx2-  (Zx)'  ' )|w  Z 


Y - ( Z-  y) 


Each  of  the  M * P numerators  require  one  add  and  two  multiplications  for  a total  of 

M * P adds  and 

2 * M * P multipl ies. 

Each  of  the  M radicals  involving  Y require  one  add,  two  multiplies,  and  one  square 
root  for  a total  of 

M adds, 

2 * M multiplies  and 
M square  roots. 

Each  of  the  M • P radicals  involving  X require  one  add,  two  multiplies,  and  one 
square  root  for  a total  of 

M ♦ P adds, 

2 ♦ M * P multiplies  and 

M * P square  roots. 
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The  total  number  of  adds  and  multiplies  for  the  five  components  is 


T^'  = M * W * (3  * P + 2) 

T,^'  = M * W * (2  * P + 1) 

The  five  components  are  combined  in  the  following  way  to  produce  the  correlation 
value  r: 


r = 


Y 


Each  of  the  M * P numerators  require  one  add  and  two  multiplications  for  a total  of 

M * P adds  and 

2 * M * P multipl ies. 

Each  of  the  M radicals  involving  Y require  one  add,  two  multiplies,  and  one  square 
root  for  a total  of 

M adds, 

2 * M multiplies  and 
M square  roots. 

Each  of  the  M ♦ P radicals  involving  X require  one  add,  two  multiplies,  and  one 
square  root  for  a total  of 

M * P adds, 

2 * M * P multiplies  and 

M * P square  roots. 
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Each  of  the  M * P denominators  require  one  multiplication  involving  the  two  radicals, 
and  finally  each  of  the  M ♦ P correlation  values  require  one  division. 

The  total  number  of  operations  involved  in  combining  the  five  components  is 


n 


M * (2  * P + 1) 
M * (5  * P + 1) 
M * P 

M * (P  + 1) 


adds 

multiplies 
divides 
square  roots. 


The  number  of  operations  used  to  calculate  the  five  components  are  combined  with 
these  results  to  produce  the  total  number  of  computer  operations  needed  to  develop 
the  M correlation  functions. 


T^  = M*|2*P+l+w*(3*P  + 2j|  adds 

T„  = M * [5  * P + 1 + W M2  * P + 1 j multiplies 


Tq  = M * P divides 

Tp  = M * (P  + 1)  square  roots. 


Computational  Overlap  Considered.  Most  of  the  computational  overlap 
occurs  in  the  calculation  of  the  five  primary  components  given  above.  The  computa- 
tional overlap  can  be  deduced  by  examining  a simple  example.  Suppose,  N = 9, 

W = 3,  and  P = 5;  then  M = 3 and  points  4,  5,  and  6 of  Y will  be  matched  to  cor- 
responding points  of  X. 

X;  Xj.  Xj. X9 

i 

i 
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The  quantities  needed  to  compute  the  five  primary  components  for  the  three  correla- 
tion functions  are  given  next. 


(Y4,  XO 


'Y  : (Y3  + Y4  + Y5) 

2,2  2 2. 

' Y : ( Y3  + Y4  + Y5) 


X : (Xj  + X2  + X3),  (X2  + X3  + XJ.  (X3  + X^  + X5),  (X^  + X5  + Xg). 


.X 


v2. 


(X5  + Xg  + X7). 

222  222  222  222 

(Xi  + X2  + X3),  (X2  + X3  + XJ,  (X3  + X^  + X5),  (X4  + X5  + Xg), 

2 2 2^ 

(X5  + Xg  + X7). 


XY:  (Xj  Y3  + X^  Y^  + X3  Yg),  (X^  Y3  + X3  Y^  + X^  Y^),  (X3  Y3  + X^  Y^  + X^  Y^), 

(X4  Y3  + Xg  Y^  + Xg  Y5),  (X5  Y3  + Xg  Y^  + X2  Yg). 

(Yj.  X5) 


V ^ • 

(Y4 

+ 

Y5 

^ Yg) 

: 

(< 

+ 

Y5 

^ Yg) 

rX  : 

ih 

+ 

X3 

+ xj. 

(X3 

+ X4 

+ Xg), 

(X4  + Xg  + Xg). 

(Xg  + 

Xe 

+ X7), 

(Xe 

+ 

X7 

+ Xg). 

iX^  : 

(X^ 

+ 

X? 

+ xj). 

(X3 

+ xj 

+ Xg). 

(X4  + X5  + Xg). 

(Xg  + 

Xe 

+ X7), 

(Xe 

+ 

A 

+ x^). 

ZXY  : 

(X2 

Y. 

+ 

X3  Ys 

+ X4 

Ye). 

(X3  Y4 

+ X4  Yg  + Xg  Yg 

).  (X4 

Y4 

+ Xg  Y 

(X5 

Y4  + 

Xe  Y5H 

X7  Yg).  (Xg  Y, 

+ X7  Y 

5^ 

Xg  Yg 
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R 


lY  : (Ys  + Ye  Yt) 

^Y^  (Ys  + ^ Y^) 

rx  : (X3  + X4  + X5),  (X4  + X5  + Xg),  (X5  + Xe  + X7),  (Xg  + X7  + Xg) 


(X7  + Xg  + Xg) 

(X3  + xj  + Xg),  (xj  + x^  + x^).  ml  + X^  + X^).  (X^  + X^  + Xg") 
222 
(X7  + Xg  + Xg) 


2XY:  (X3  Yg  + X4  Ye  + Xg  Y7).  (X^  Y^  + X^  Y^  + X^  Y^).  (X^  Y^  + Xg  Y^  + X^  Y. 

(Xg  Yg  + X,  Yg  + Xg  Y^).  (X,  Yg  + Xg  Yg  + Xg  Y,) 


The  2 Y value  for  the  first  correlation  function  can  be  developed  by  W adds, 
and  the  remaining  ZY  values  for  the  (M-1)  other  correlation  functions  can  be  de- 
veloped by  two  adds  apiece  for  a total  of  W + 2 * (M-1)  adds.  The  ZY  value  for 
R(Ys  ,X  5 ) is  computed  by  subtracting  Y3  from  and  adding  Y6  to  the  ZY  value  for 
R(Y4,X4).  The  ZY  value  for  R(Y6,X6)  is  computed  by  subtracting  Y4  from  and 
adding  Y?  to  the  Z Y value  for  RfYs  ,Xs ).  The  Z Y^  values  are  developed  in  the  same 
manner;  the  results  for  the  Z Y and  ZY  components  are 

^jY:  W + 2 * (M-1)  adds 

; W + 2 * (M-1)  adds 

W + 2 * (M-1 ) multi  pi ies 

The  ZX  values  can  be  computed  even  more  efficiently  than  the  ZY  values. 
Cbnsider  the  P = five  ZX  values  for  R(Y4,X4).  The  first  can  be  developed  with  W 
adds,  and  the  remaining  (P-1)  values  can  be  developed  with  two  adds  apiece  for  a 
total  of  W + 2*  (P-1)  adds.  Note  that  (P-1)  of  the  ZX  values  for  RfYs , Xs ) are 
identical  to  those  of  R(Y4  ,X4 ),  and  the  remaining  ZX  value  can  be  developed  by  two 
adds.  The  total  number  of  adds  for  ZX  is  W -H  2 * (P-1)  + 2*  (M-1).  The  ZX’  values 
a-e  developed  in  the  same  way;  the  results  for  the  ZX  and  ZY’  values  are 


W + 2 * (P-1)  + 2 * (M-1)  adds 


W + 2 * (P-1)  + 2 * (M-1)  adds 


W + 2 * (P-1)  + 2 * (M-1)  multiplies 


There  is  no  computational  overlap  for  the  SXY  components  within  a specific 
correlation  function,  but  there  is  computational  overlap  between  correlation  func- 
tions. The  first  correlation  function  will  require  W * P adds  and  W * P multiplies  to 
compute  the  P 2XY  values.  Each  of  the  succeeding  (M-1)  correlation  functions  will 
require  2 * P adds  and  P multiplications  for  a total  of 

XxY.  w*p  + 2*P*  (M-1)  adds 

W * P + P * (M-1)  multiplies 


The  total  number  of  required  adds  and  multiplies  for  the  five  components  are 

TA  = W * (P+4)  + 2 * (M-1)  * (P+4)  + 4 * ( P-1) 

TM  = W * (P+2)  + (M-1)  * (P+4)  + 2 * (P-1) 

Except  for  the  radical  involving  X,  the  number  of  individual  calculations 
(where  the  five  components  are  combined  to  produce  correlation  values)  is  identical 
to  the  number  of  operations  for  the  radical  involving  X is 

M + P - 1 adds 
2 * M + 2 * (P-1)  multiplies  and 

M + P - 1 square  roots. 


The  total  number  of  computer  operations  involved  in  combining  the  five  components 
is 
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it 


T;\  = M * (P+2)  + (P-l) 
Tm=M*(3*P+4)+2*  (P-1 ) 

II 

Tq  = M * P 

Tr"  = 2 * M + (P-l ) 


The  numbers  of  operations  used  to  calculate  the  five  components  are  com- 
bined with  the  results  above  to  produce  the  total  number  of  computer  operations 
needed  to  develop  the  M correlation  functions. 


= [W  + 3 * (M-bljJ  * P + (4  * W + 10  * M - 13) 

T|^  = ['W+4*M  + 3j*P  + 2*  [w+8*(M-l)J 

Tp  = M * P 

Tp  = 2 * M + (P  - 1) 

The  values  R^,  R|^,  Rq,  and  Rj^  are  the  multiplicative  increases  in  compute 
operations  when  computational  overlap  is  ignored.  The  value  R^  (increase  in  number 
of  adds)  is  computed  by  calculating  the  ratio  of  the  two  T^^  values;  the  other  ratios 
are  calculated  in  the  same  manner. 

Ra  (3  * P 2)  * W 

(3  * P + 10)  (P  + 4) 

P-M  (2  * P + 1)  * W 

N 

% =’ 

Rr  •>.  P/j 

if,  for  example,  N is  very  large  compared  to  the  window  size  W and  if  P is 
also  large,  then  R^^-vW,  Rm~W/2  and  R'''P/2. 
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